Module 6: Random Effects & Linear Mixed Models

Linear Mixed Models

Example 6.3: Knee Brace

Researchers study 4 knee braces: Control, Air Armor 2, Breg, Don Joy

10 high school running backs with previous ACL injury participate. Each player: Wears all 4 braces, Order randomized

Time to complete 40-yard agility course recorded

Example 6.3: Skeleton ANOVA

Do we care only about these 10 running backs?

Source of Variation	DF

When one factor is fixed and one is random

Linear Mixed Models (LMM)

Mixed models come about when one factor is fixed and one factor is random. A mixed model should be used to study the effects of one factor on the population mean and the effects of another factor on the population variance.

Statistical Effecs Model

\[y_{ij} = \mu + \tau_i + p_j + \epsilon_{ij} \text{ with } p_j \sim \text{iid } N(0,\sigma_{blk}^2) \text{ and } \epsilon_{ij} \sim \text{iid } N(0,\sigma_\epsilon^2)\]

for \(i=1,2,3; j=1,2,…,10\)

where:

\(y_{ij}\) is the time to complete a 40-yard agility course for the \(j^{th}\) running back wearing the \(i^{th}\) knee brace
\(\mu\) is the overall mean time to complete the agility course
\(\tau_i\) is the fixed effect of the \(i^{th}\) knee brace
\(p_j\) is the random effect of the \(j^{th}\) running back (block)
\(\epsilon_{ij}\) is the random error term associated with the \(j^{th}\) running back wearing the i^th knee brace

Scope of Inference

Note

If blocks were fixed: Inference limited to these 10 running backs
If blocks are random: Inference extends to all similar running backs

The Data

library(tidyverse)
kneebrace_data <- read_csv("data/06-kneebrace-data.csv") |> 
  mutate(across(RunningBack:KneeBrace, as.factor))
head(kneebrace_data, n = 12)

# A tibble: 12 × 3
   RunningBack KneeBrace AgilityTime
   <fct>       <fct>           <dbl>
 1 1           Control          19.1
 2 1           AirArmor         18  
 3 1           Breg             18.0
 4 1           DonJoy           17.9
 5 2           Control          19  
 6 2           AirArmor         18.4
 7 2           Breg             18.2
 8 2           DonJoy           17.7
 9 3           Control          19.9
10 3           AirArmor         18.7
11 3           Breg             18.6
12 3           DonJoy           18.7

R: Fitting a LMM

library(lme4)
kneebrace_mod <- lmer(AgilityTime ~ KneeBrace + (1 | RunningBack), data = kneebrace_data)
summary(kneebrace_mod)

Linear mixed model fit by REML ['lmerMod']
Formula: AgilityTime ~ KneeBrace + (1 | RunningBack)
   Data: kneebrace_data

REML criterion at convergence: 1.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.72607 -0.47036 -0.03709  0.53329  1.71255 

Random effects:
 Groups      Name        Variance Std.Dev.
 RunningBack (Intercept) 0.15453  0.3931  
 Residual                0.01965  0.1402  
Number of obs: 40, groups:  RunningBack, 10

Fixed effects:
                 Estimate Std. Error t value
(Intercept)      18.24200    0.13198  138.22
KneeBraceBreg    -0.01000    0.06269   -0.16
KneeBraceControl  0.97500    0.06269   15.55
KneeBraceDonJoy  -0.16800    0.06269   -2.68

Correlation of Fixed Effects:
            (Intr) KnBrcB KnBrcC
KneeBracBrg -0.237              
KneBrcCntrl -0.237  0.500       
KneeBrcDnJy -0.237  0.500  0.500

JMP: Fitting a LMM

Analyze > Fit Model

Testing Hypotheses About the Treatment Means

\[H_0: \tau_1=\tau_2=\tau_3=\tau_4 \text{ vs } H_A: \text{ at least one } \tau_i \text{ differs}\]

library(lmerTest)
anova(kneebrace_mod)

Analysis of Variance Table
          npar Sum Sq Mean Sq F value
KneeBrace    3 8.2015  2.7338  139.14

Estimating Treatment Means (balanced design)

Recall \(y_{ij} = \mu + \tau_i + p_j + \epsilon_{ij}\)

In the mixed model:

\(E[y_{ij}] = \hat\mu+\hat\tau_i\)
\(Var(y_{ij}) = \hat\sigma_{blk} + \hat\sigma_\epsilon^2\).

Thus, the estimated mean and variance are given by:

Mean: \(\hat\mu_i=\hat\mu+\hat\tau_i\)
SE of Mean: \(SE(\hat\mu_i)=\sqrt{\frac{\hat\sigma_{blk}^2+\hat\sigma_\epsilon^2}{r}}\)
SE of Diff in Means: \(SE(\hat\mu_i-\hat\mu_{i'})=\sqrt{\frac{2}{r}\sigma_\epsilon^2}\)

Estimating Treatment Means

library(emmeans)
library(multcomp)
# emmip(kneebrace_mod, ~ KneeBrace, CI = T)
emmeans(kneebrace_mod, specs = ~KneeBrace) |> 
  cld(Letters = LETTERS, decreasing = T, adjust = "tukey")

 KneeBrace emmean    SE   df lower.CL upper.CL .group
 Control     19.2 0.132 10.7     18.8     19.6  A    
 AirArmor    18.2 0.132 10.7     17.8     18.6   B   
 Breg        18.2 0.132 10.7     17.8     18.6   B   
 DonJoy      18.1 0.132 10.7     17.7     18.5   B   

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 
Conf-level adjustment: sidak method for 4 estimates 
P value adjustment: tukey method for comparing a family of 4 estimates 
significance level used: alpha = 0.05 
NOTE: If two or more means share the same grouping symbol,
      then we cannot show them to be different.
      But we also did not show them to be the same.

Effect on Standard Errors

Treating blocks as random (instead of fixed) (1) Increases SE of treatment means (2) Does NOT change SE of treatment differences

Running back Fixed (LM)

kneebrace_lm <- lm(AgilityTime ~ KneeBrace + RunningBack, 
                      data = kneebrace_data)
emmeans(kneebrace_lm, specs = ~KneeBrace)

 KneeBrace emmean     SE df lower.CL upper.CL
 AirArmor   18.24 0.0443 27    18.15    18.33
 Breg       18.23 0.0443 27    18.14    18.32
 Control    19.22 0.0443 27    19.13    19.31
 DonJoy     18.07 0.0443 27    17.98    18.16

Results are averaged over the levels of: RunningBack 
Confidence level used: 0.95

emmeans(kneebrace_lm, specs = ~KneeBrace) |> pairs()

 contrast           estimate     SE df t.ratio p.value
 AirArmor - Breg       0.010 0.0627 27   0.160  0.9985
 AirArmor - Control   -0.975 0.0627 27 -15.553 <0.0001
 AirArmor - DonJoy     0.168 0.0627 27   2.680  0.0565
 Breg - Control       -0.985 0.0627 27 -15.713 <0.0001
 Breg - DonJoy         0.158 0.0627 27   2.520  0.0792
 Control - DonJoy      1.143 0.0627 27  18.233 <0.0001

Results are averaged over the levels of: RunningBack 
P value adjustment: tukey method for comparing a family of 4 estimates

Running back Random (LMM)

kneebrace_lmer <- lmer(AgilityTime ~ KneeBrace + (1 | RunningBack), 
                      data = kneebrace_data)
emmeans(kneebrace_lmer, specs = ~KneeBrace)

 KneeBrace emmean    SE   df lower.CL upper.CL
 AirArmor    18.2 0.132 10.7     18.0     18.5
 Breg        18.2 0.132 10.7     17.9     18.5
 Control     19.2 0.132 10.7     18.9     19.5
 DonJoy      18.1 0.132 10.7     17.8     18.4

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95

emmeans(kneebrace_lmer, specs = ~KneeBrace) |> pairs()

 contrast           estimate     SE df t.ratio p.value
 AirArmor - Breg       0.010 0.0627 27   0.160  0.9985
 AirArmor - Control   -0.975 0.0627 27 -15.553 <0.0001
 AirArmor - DonJoy     0.168 0.0627 27   2.680  0.0565
 Breg - Control       -0.985 0.0627 27 -15.713 <0.0001
 Breg - DonJoy         0.158 0.0627 27   2.520  0.0792
 Control - DonJoy      1.143 0.0627 27  18.233 <0.0001

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 4 estimates